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ABSTRACT 


In  the  computation  of  discontinuous  solutions  of  hyperbolic  systems  of  conservation  laws, 
the  recently  developed  ENO  (Essentially  Non-Oscillatory)  schemes  appear  to  be  very  useful. 
However,  they  are  computationally  costly  compared  to  simple  central  difference  methods.  In 
this  paper  we  develop  a  filtering  method  which  uses  simple  central  differencing  of  arbitrarily 
high  order  accuracy,  except  when  a  novel  local  test  indicates  the  development  of  spurious 
oscillations.  At  these  points,  generally  few  in  number,  we  use  the  full  ENO  apparatus, 
maintaining  the  high  order  of  accuracy,  but  removing  spurious  oscillations.  Numerical  results 
indicate  the  success  of  the  method.  We  obtain  high  order  of  accuracy  in  regions  of  smooth 
flow  without  spurious  oscillations  for  a  wide  range  of  problems  and  a  significant  speed  up  of 
generally  a  factor  of  almost  three  over  the  full  ENO  method. 
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1  Introduction 


Recently,  a  new  class  of  high  order  schemes  related  to  numerical  calculations  of  linear  and 
nonlinear  systems  of  conservation  laws  has  been  developed.  This  new  class  of  methods 
uses  an  explicit  TVD  Runge-Kutta  Multistage  time  discretization  together  with  high  order 
Essentially  Non  Oscillatory  (ENO)  spatial  discretization  methods.  Details  of  these  methods 
can  be  found  in  [6,  7]  and  in  references  quoted  therein.  These  methods  are  applied  to  solve 
numerically  hyperbolic  systems  of  conservation  laws 

ut  +  (f(u))x  =  0,  (1) 

u(x,0)  =  u0(x),  (2) 

to  be  solved  for  t  >  0  and  x  in  some  interval  Q  with  appropriate  boundary  conditions.  An 
example  of  such  a  system  of  conservation  laws  ;s  given  by  the  Euler  equations  of  compressible 
gas  dynamics  for  which 

f(u)  =  uu  +  (0,p,t;p)T  (3) 

and  u  =  ( p ,  q,p ),  p  is  density,  q  is  momentum,  v  is  velocity,  and  p  is  the  pressure.  A  similar 
example  for  two  dimensional  flows  will  be  considered.  The  two  dimensional  version  of  the 
equation  (1)  now  has  different  fluxes  for  each  space  dimension.  We  have: 

u«  +  (f(«)).  +  (g(u)v  =  0,  (4) 

u(x,3/,0)  =  u0(x,y),  (5) 

to  be  solved  for  t  >  0,  (x,y)  6  D,  some  compact  set,  with  appropriate  boundary  conditions. 
The  fluxes  are  f(u)  =  uxu-|-(0,p,  0,u*p)T  and  g(u)  =  uvu  +  (0,0 ,p,vvq)T,  respectively,  where 
u  =  (p>  9xi  9v>  e)T •  numerical  experiments,  we  approximate  the  solution  of  equations  (  1) 
or  (  4)  by  using  point  values.  That  is,  u (xj,tn)  is  approximated  by  u",  given  a  regular 
triangulation  of  the  domain  fi.  In  this  paper,  only  a  line  by  line  discretization  will  be 
considered,  restricting  the  shape  of  domain  fl  to  regions  which  can  be  mapped  onto  squares 
or  rectangles. 

The  TVD  time  discretization  performed  is  the  one  introduced  in  [6,  7].  The  method  is 
explicit  and  relatively  easy  to  program.  Such  algorithms  can  be  briefly  described  as  follows: 

un+i/m  =  g[aau»+V~  +  /3iifcAiL(u"+fc/"*)],  (6) 

k=0 

where  m  is  the  number  of  stages  to  move  the  solution  from  time  t  to  t  +  At.  Generally, 
second  to  sixth  order  methods  are  investigated.  The  coefficients  and  (3iilc  are  calculated 
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so  as  to  improve  the  CFL  number  in  order  to  minimize  the  number  of  time  iterations.  In 
the  class  of  method  defined  by  (  6),  the  CFL  coefficient  A  must  satisfy: 


where  A0  is  the  maximum  allowable  value  for  the  forward  Euler  method  (see  [6]).  In 
particular,  it  is  possible  to  obtain  a  third  order  accurate  TVD  time  discretization  method 
with  a  CFL  number  of  1.  This  is  only  slightly  reduced  (to  |)  for  the  fourth  order  method. 
In  addition,  it  is  possible  to  derive  a  class  of  time  discretizations  that  require  the  evaluation 
of  L(un+(‘-1A/,fc)  onV  so  that  /?o,  -  o.  This  process  reduces  the  otuiage  xequiremcut 

significantly.  However,  such  procedure  is  possible  only  for  methods  of  up  to  third  order 
accuracy.  For  higher  order  methods,  several  evaluations  of  the  operator  L  are  needed  to 
enforce  the  TVD  property;  see  [6,  7]  for  more  details. 

For  the  space  discretization,  we  use  high  order  ENO  methods  to  approximate  L(U)  using 
conservation  form. 


-L{uj)  = 


f(UJ— r+l/2j  •••juj+»+l/2)  -  f(uj-r-l/2>"-)Uj+»-l/2) 


For  multidimensional  operator,  —  L  =  f*  +  gv,  —  f*  is  performed  as  in  (  7),  and  — g„  is 
approximated  analogously.  For  systems  of  equations,  a  field  by  field  decomposition  is  used. 
We  calculate  at  each  point  the  eigen  decomposition  of  different  fluxes,  evaluating  the  rth 
order  accurate  interpolating  polynomial  that  approximates  the  fluxes  in  each  field,  and  then 
recover  each  vector  field  in  the  physical  space  by  inverse  decomposition.  In  many  cases  (Euler 
equations  for  example),  the  decomposition  in  each  field  uses  the  left  eigenvectors  L*,  so  that 
ffc,j+i/2  =  Lfc.A(uJJ+i),  where  A  is  the  Roe  matrix  for  Vf(u)  (see  [10]  for  more  details).  The 
ENO  algorithm  is  based  on  a  Newton  interpolating  polynomial  using  an  adaptative  stencil. 
That  is,  instead  of  considering  a  polynomial  interpolant  using  a  fixed  centered  or  a  fixed 
upwind  stencil,  we  derive  an  interpolating  polynomial  minimizing  the  successive  undivided 
differences.  This  process  limits  oscillations,  thus  the  name  of  the  method.  In  the  case  of  a 
shock  discontinuity,  this  method  works  quite  well,  leading  to  sharp  transitions  over  a  few 
points.  However,  smearing  of  linear  (e.g  contact)  discontinuities  may  occur.  In  this  case,  a 
particular  treatment  using  subcell  resolution  or  artificial  compression  is  available  to  sharply 
resolve  large  transitions.  The  interested  reader  can  find  more  details  in  [4,  7,  12],  we  do  not 
use  this  improvement  here. 

The  reader  has  probably  already  realized  that  such  methods  require  a  considerable 
amount  of  computational  time  when  multidimensional  systems  are  investigated.  It  would  be 
interesting  to  use  high  order  methods  derived  from  simple  spatial  discretization  most  of  the 
time,  and  use  ENO  methods  only  when  spurious  oscillations  appear.  Of  course,  it  is  well 
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known  that  such  oscillations  will  generally  occur  for  nonlinear  equations  even  with  smooth 
initial  conditions. 

In  order  to  deal  with  problems  in  which  singularities  occur  but  still  use  simple  interpola¬ 
tion  techniques,  hence  simple  finite  difference  methods,  we  introduce  a  postprocessing  step 
that  detects  existing  singularities  and  spurious  oscillations  and  corrects  them  if  necessary. 
Such  a  method  based  on  a  postprocessing  step  has  been  introduced  in  [1,  11].  It  relies  on  a 
simple  type  of  correction  by  pushing  points  up  or  down  up  to  an  acceptable  level  so  that  the 
global  solution  satisfies  both  total  variation  diminishing  (TVD)  and  conservation  properties. 
By  modifying  this  type  of  filter,  the  third  author  in  [1]  was  able  to  prove  convergence,  in 
[11],  to  the  physical  solution  for  one  dimensional  nonlinear  conservation  laws.  His  proof 
relies  on  compensated  compactness  arguments  using  Young  measures.  Some  references  can 
be  found  in  [11]  and  in  other  papers  quoted  therein.  More  importantly,  an  extremely  simple 
but  useful  TVD  algorithm  which  works  very  well  in  practice  was  developed  in  [1,  11]. 

In  this  paper,  we  define  a  new  class  of  filters  enforcing  uniformly  high  order  of  accuracy 
without  allowing  significant  spurious  oscillations.  Section  2  is  devoted  to  a  detailed  presen¬ 
tation  of  our  method  .  Section  3  offers  several  numerical  examples  in  one  and  two  space 
dimensions  for  both  scalars  and  systems  of  conservation  laws.  Section  4  will  conclude  this 
paper  and  propose  several  different  approaches  for  solving  nonlinear  problems  in  general 
domains  using  finite  element  grids. 

2  High  order  uniform  filtering  methods 

We  now  outline  our  filtering  method.  We  first  only  consider  spatially  centered  differences. 
This  process  leads  to  a  simple  scheme  with  low  computational  cost  for  the  evaluation  of 
the  numerical  fluxes.  Then,  from  the  solution  which  has  been  evaluated  from  this  basic 
simple  scheme,  we  perform  another  step  that  filters  the  numerical  oscillations  by  using  high 
order  accurate  ENO  interpolation.  We  then  use  a  high  order  TVD-RK  time  scheme.  The 
algorithm  is  thus  simply: 

•  For  i  =  l,m  do 

1.  Approximate  the  equation  (  1)  or  (  4)  using  the  basic  scheme 

vn+i/m  _  £)(„«  .  .  .  >  un+(i-l)/m^  (g) 

where  D  is  the  numerical  operator  that  performs  the  m- multistage  TVD-RK 
algorithm  as  in  [6]  together  with  centered  spatial  differences  approximating  the 
operator  L  in  (  7)  in  conservation  form. 
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2.  If  vn+‘/m  has  spurious  oscillations  Then  we  correct  it  using  the  filter: 

un+»/m  _  p^unj  .  .  .  }  un+(*-2)/m^  vn+(«-l)/m^  (g) 

where  F  is  the  numerical  operator  that  uses  the  same  time  discretization  algorithm 
as  the  basic  scheme  together  with  a  high  uniform  ENO  type  filter  in  evaluating 
numerical  fluxes. 

•  End  For. 

We  define  high  order  2 pth  approximations  as  follows: 

•  Second  Order:  ( p  =  1) 


fj+l/2  - 


f(u,-+i)  +  f(u  j) 


•  Fourth  Order:  ( p  =  2) 


fj+1/2  =  ^(f(ui+0  +  ^(uj))  ~  y2^Uj+2  +  f(u3-1))J 


(10) 


(11) 


•  Higher  Order  Vp:  It  is  a  simple  matter,  using  Richardson’s  extrapolation  to  construct 
and  obtain  arbitrary  high  order  accurate  centered  difference  methods.  This  has  been 
done  in  many  places,  e.g  [5].  We  will  have: 

p 

fj+1/2  =  53  A^(ui+»)> 

»=-p+i 

where  ^_<+i  =  Pi  for  i  =  —  p  +  1,  ...,p  —  1, 
and  ]T  &  =  1- 

♦=-p+i 

For  the  time  discretization,  we  consider  the  TVD  Runge-Kutta  type  idea  introduced  in 
[6,  7],  that  is 

•  For  second  order  method:  (p  =  1) 

un+l/2  =  u"  +  AtL(U"), 

ur 


,n+l  _  Iun+l/2  4.  I„*»  _l  A#r('iin+1/2 
2 


+  lun  +  AtL(un+1/2). 

ml 


(12) 

(13) 


•  For  Third  order  method:  (p  =  3/2) 

un+l/3  _  un  +  AtZ(Un), 

un+2/3  =  +  iUn+l/3  +  7AtX(Un+l/3), 

4  4  4 

Un+1  _  Iu"  +  Hu"+2/3  +  -/\tL(un+2/3). 

3  3  3 
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(14) 

(15) 

(16) 


•  Higher  order  methods  of  this  type  are  described  in  [6]  up  to  sixth  order. 


The  filter  step  changes  the  centered  differences  spatial  approximation  to  the  more  stable 
ENO  approximation  of  fluxes.  As  building  block,  we  use  either  the  Roe  scheme  (see  [10,  7]) 
which  admits  expansion  shocks  at  sonic  points  or  the  Local  Lax-Friedrichs  decomposition 
of  the  fluxes  at  such  points,  see  [7].  Both  building  blocks  first  decide  on  the  initial  stencil 
that  respects  the  local  characteristic  direction  and  then  evaluate  the  polynomial  interpolant 
using  an  adaptative  stencil.  This  stencil  is  chosen  in  order  to  minimize  derivatives  of  the 
interpolating  polynomial.  The  algorithm  for  computing  the  numerical  fluxes  in  the  filtering 
step  is  precisely  the  algorithm  2.3,  in  [7].  Furthermore,  in  order  to  still  get  a  globally 
conservative  scheme,  backward  and  forward  corrections  are  performed.  This  lead  to  these 
four  possible  approximations  of  u"+*^m  after  the  filtering  step: 


n+i/m 

u; 

n+t/m 

u; 

un+t/m 

u"+,/m 


= 

=  *(uj, 

= 

= 


-.ll;+<i''l)/”)+^(fst/2-fn/2)l 


for  i  =  1  ,...,m.  The  operator  A  represents  some  linear  combination  of  un+fc/m,  for  k  — 
0,  ...,i  —  1,  given  for  example  by  the  equations  (  12), (  13),  or  (  14), (  15), (  16),  and  fc,  imo 
are  the  computed  fluxes  using  centered  differences  or  ENO  interpolants,  respectively.  If  the 
evaluation  of  f/Hi/2  needed,  then  a  correction  is  performed  on  m ,  where  j  is  a  multiple 
index  in  the  case  of  several  dimensions.  Hopefully,  for  regular  grids,  backward  corrections 
will  only  occur  at  the  first  column  and  first  row  of  the  computational  domain  fi.  Hence,  at 
interior  grid  points,  only  forward  corrections  of  the  fluxes  are  performed.  For  example,  in 
two  dimensions,  an  initialization  step  in  the  filtering  step  is  done  for  the  first  column  and 
first  row,  say  i  =  i0,  and  j  =  jo,  leading  to  corrections  of  the  fluxes  and  £+1/2, j  for 

i  >  i0  and  j  >  j0  only. 

In  the  multidimensional  case,  this  is  done  separately  in  the  x  and  y  directions.  The 
global  scheme  is  therefore  as  simple  as  for  the  one  dimensional  case.  For  systems  of  conser¬ 
vation  laws,  the  centered  approximation  is  performed  for  each  component  of  the  fluxes.  The 
ENO  interpolant,  when  needed,  must  be  evaluated  in  each  characteristic  field.  To  do  this, 
we  follow  the  algorithm  4.1  in  [7].  That  is,  we  evaluate  the  gradient  of  each  fluxes  using 
Roe  averages  for  the  unknown  A j+i/3  =  §£^=11*°^,!  where  =  R(uJ+i,  u j)  is  the  Roe 

average  of  Uj  and  uJ+i  (see  [10]).  Denoting  the  left  and  right  eigenvectors  of  Aj+1/2  by 


5 


Ly-i-i/2,  R£+i/2i  P  =  1,  ...,n,  where  n  is  the  number  of  equations,  we  correct  each  character¬ 
istic  flux  using  the  divided  differences  projected  in  each  field  B =  L^jy2.Bj+1/2, 

where  the  Bj+i/2  are  the  successive  divided  differences  of  fj+1/2  needed  in  the  evalua¬ 
tion  of  the  interpolating  polynomial.  We  then  reconstruct  each  flux  in  the  physical  space 
fj+1/2  =  £p=i  fj+i/jB^+jyj.  At  sonic  points,  we  replace  the  Roe  building  block  by  the  LLF 
decomposition  of  the  fluxes.  Also,  instead  of  natural  divided  differences  of  the  fluxes,  we  ap¬ 
ply  the  formula  (2.11a)  and  (2.11b)  of  [7]  in  each  field  * -iung  ^i/2uj-i/2)> 

and  fj^i/2  =  |(fj-i/2  +  where  A =  m^n  (uj-i> u^)  case  °f  convex 

fluxes  or  genuinely  nonlinear  fields.  For  more  details,  see  [7]. 

The  most  interesting  and  important  part  of  this  method  consists  in  the  way  large  tran¬ 
sition  areas  are  detected.  For  the  basic  idea,  we  use  concepts  related  to  a  front  capturing 
method  which  has  been  introduced  in  [8]  in  the  case  of  combustion  type  problems  dealing 
with  Hamilton-Jacobi  type  equations  and  in  many  other  applications  where  fronts  must  be 
located  with  high  accuracy.  For  conservation  laws,  we  will  say  that  from  time  t  to  time  t  +  At 
the  solution  has  changed  considerably,  if  the  normal  at  this  point  to  the  solution  surface  at 
these  two  different  times  has  rotated  by  an  angle  exceeding  some  preset  value  a,  |a|  <  |^|. 
In  numerical  computations,  spatial  derivatives  of  are  evaluated  using  backward  or 

forward  derivatives,  so  that  the  change  of  the  normal  is  described  via  this  formula: 


A±u"A±u”+1  < 


(17) 


where  the  ±  tests  are  performed  to  detect  fast  transition  locations  for  Uj+1/2  or  u^/a, 
respectively,  and  A±u”+1  =  ±(uj±j  —  Uj)  are  for  the  forward  and  backward  differences. 

If  this  inequality  is  satisfied  with  the  value  of  the  parameter  —j<a<j  then  the 
correction  step  is  performed.  The  parameter  a  is  introduced  to  restrict  more  general  changes 
of  the  normal  to  the  surface  from  different  times.  In  the  numerical  examples  introduced  in 
next  section,  it  has  been  necessary  to  tune  this  coefficient  in  order  to  obtain  the  smallest 
possible  number  of  corrections  without  introducing  oscillations  in  the  numerical  solution. 
If,  instead,  we  let  cos  a  =  0,  more  corrections  must  be  performed.  Basically,  we  also  have 
to  correct  when  an  extremum  point  already  exists  at  time  t.  This  test  can  be  avoided 
when  larger  values  of  the  parameter  cos  a  are  used.  This  process  has  greatly  reduced  the 
number  of  corrections  when  two  dimensional  Euler  equations  for  compressible  gas  dynamics 
with  shock-turbulence  interaction  has  been  studied.  However,  it  has  been  difficult  to  adjust 
this  parameter  to  optimize  the  global  algorithm  (computational  time).  Nevertheless,  less 
than  40%  of  the  grid  points  required  corrections  for  cos  a  =  0.1.  Further  optimization  will 
probably  lead  to  even  fewer  corrections. 
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The  value  cos  a  =  1  implies  that  the  filtering  scheme  is  applied  everywhere,  while,  cos  a  = 
—  1  leads  only  to  the  simple  basic  centered  difference  scheme.  The  particular  value  of  the 
parameter  cos  a,  for  which  the  right  hand  side  of  (  17)  is  zero,  implies  that  a  change  in  sign 
of  a  partial  derivative  occured  from  time  t  to  t  -f  At. 

We  must  add  that  no  significant  improvement  of  the  solution  is  obtained  when  the  value 
of  cos  a  approaches  1.  This  is  probably  due  to  the  fact  that  the  numerical  error  is  very  small 
in  regions  where  the  physical  solution  is  smooth  for  which  centered  differences  are  used. 
Furthermore,  as  shown  in  the  numerical  examples  below,  there  is  no  visible  propagation  of 
the  numerical  error  due  to  the  filtering  step. 

To  conclude  this  section,  we  write  down  clearly  the  algorithm  for  multidimensional  sys¬ 
tems  of  conservation  laws: 

•  For  i  =  l,  ...,m  Do 

—  For  Ij  =  1,  Do,  j  =  1,  ...,nd,  nd  is  the  spatial  dimension, 

*  Compute  j  using  the  basic  centered  difference  scheme: 

n+i/m  _  A(..n  n+(t-l)/mx 

Uh . /nd  “  ) 

At 

^  Ax  I  . /„d  fyt  -1/2,/a 

+  —  +  . /nd+1/2  “  f/l /nd-l/z) 

-  End  For. 

-  For  Ij  =  1, ..., Nj  Do  ,  j  =  1, ... ,nd : 

*  Compute  the  normal  to  the  surface  at  Lime  and  and  test 

whether  the  directional  change  of  the  normals  exceeds  the  angle  a: 

Axh  „"+(»-*)/"»  A,xh  n+i/m 

U/l . Ini  a±  U/l . /nd  , 


A  *'nd  ,  A  *'nd 

U/l . /nd  U/l . /nd 


<  1  +  cos  a 


AiL 


ax/i  1In+(*-1)/m 

1+(  ,  )>  + ■■■  +  (-*  ^  f 

1  + +  -  +  ('  *A  I: . '-')*) 


*  If  these  tests  are  satisfied  Then  compute  the  ENO  interpolant  in  each  field 
and  correct  the  solution. 

-  End  For. 

•  End  For. 

The  number  of  tests  to  be  performed  depends  on  the  dimension  of  the  system  (=  ns)  and  of 
space  dimension  (=  nd ).  In  general,  we  will  have  ns  *nd  tests  to  be  performed  so  that  each 
component  of  the  fluxes  in  each  direction  can  be  modified.  Therefore  this  postprocessing 
may  imply  a  high  computational  cost  if  the  code  is  not  well  implemented. 

3  Numerical  Results 

We  now  apply  our  algorithm  to  several  test  problems  dealing  with  linear  and  nonlinear 
hyperbolic  systems  of  conservation  laws.  We  will  focus  our  attention  in  determining  precisely 
the  order  of  accuracy  of  our  filtering  method  and  studying  the  propagation  of  the  local  error. 
We  will  also  visualize  the  numerical  solution  when  linear  (contact)  or  nonlinear  (shock) 
discontinuities  appear  and  compare  the  computational  time  of  the  filtering  method  versus 
the  associated  unfiltered  ENO  scheme. 


3.1  Example  I: 

As  first  test  problem,  we  want  to  verify  that  our  method  is  uniformly  high  order  accurate 
even  at  extremum  points.  To  do  so,  we  consider  the  simple  linear  equation 


with  initial  condition 


ut  +  ux  =  0, 


u(x,0)  =  cos  2IIx, 

to  be  solved  for  t  >  0  and  x  £  [0, 1]  with  periodic  boundary  conditions  at  0  and  1.  Numer¬ 
ically,  we  discretize  the  interval  [0,1]  and  let  x^  =  i Ax,  for  i  =  0,...,n.  The  exact  solution 
u(xj,  tn)  is  then  approximated  by  pointwise  values  u”  for  which  we  set  =  u(xi,  0)  from  the 
initial  condition.  The  exact  solution  is  calculated  by  the  method  of  characteristics  so  that  the 
local  numerical  error  and  order  of  accuracy  can  be  estimated.  In  the  numerical  experiments, 
we  fixed  n  =  40  and  ran  the  program  for  one  period  of  time,  e.g  t  =  1.  The  number  of  cor¬ 
rections  per  each  time  substep  was  never  higher  than  4.  Moreover,  these  corrections  occur  at 
extremum  points.  Table  1  shows  the  local  error  at  points  x^j  —  4;. Ax.  Table  2  describes  the 
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global  order  of  accuracy  in  L 1  and  L°°  norms  for  the 
x-location  (3-2)FM  Local  Error  (3-4)FM 
~0.  7.15  *  10"2 

0.1  2.60  *  10“a 

0.2  1.96*  10“2 

0.3  3.27  *  10"2 

0.4  1.25  *  10-3 

(FM)'  T3 - T24.10-> - 

0.6  2.71  *  10"2 

0.7  2.02  *  10-2 

0.8  3.31  *  10~3 

0.9  1.31  *  10-s 

1.  6.29  *  10~J 


(3-2),  (3-4),  and  (3-6)  filtering  methods 


Local  Error 

(3-6)FM  Local  Error 

5.69  *  10-4 

1.34  *  10"4 

1.70*  10~4 

1.01  *  10"4 

1.26*  10~4 

3.36  *  10"5 

6.91  *  10~4 

4.75  *  10-5 

1.05  *  10~4 

1.01  *  lO"4 

5.59  *  10~4 

1.34*  10“4 

1.69*  10~4 

1.01  *10"4 

1.24*  10~4 

3.36  *  lO"5 

7.32  *  10“5 

4.74*  10"5 

1.10*  10~4 

1.27*  lO'4 

5.58  *  10~4 

1.32  *  10^ 

Scheme 

LUnorm 

L°°-norm 

(3-2)FM 

2.08 

1.75 

(3-4)FM 

4.11 

3.46 

(3-6)FM 

3.46 

3.26 

!able  1. 


Table  2. 


These  results  are  indeed  in  agreement  with  what  we  should  expect.  Uniform  high  order 
is  preserved  even  at  corrected  points. 


3.2  Example  II: 

We  now  extend  Example  I  to  two  dimensions.  As  described  in  the  algorithm,  a  dimension 
by  dimension  approach  is  used  to  solve  the  linear  equation 

V-t  +  Hi  +  Uy  =  0, 

u(x,y,0)  =  cos  2II(x  +  y)  and  0  <  x,y  <  1, 

and  again,  periodic  boundary  conditions  in  both  x  and  y  variables  are  assumed.  We  discretize 
the  square  domain  fi  =  [0, 1]  x  [0, 1]  and  denote  by  the  vertices  of  coordinates  x<  =  iAx, 
and  yj  =  jAy,  for  i  =  0,  ...,n  and  j  =  0,  ...,m.  We  choose  n/mso  that  the  problem  is 
really  two  dimensional.  The  exact  solution  of  this  equation  can  be  easily  calculated  using 
the  change  of  variables  £  =  x  +  y  leading  to  a  one  dimensional  problem.  Hence,  the  exact 
solution  is  merely 

u(x,  y,  t )  =  cos  (2II(x  +  y  -  t)). 

Numerically,  we  use  m  =  40,  and  n  =  30  so  that  Ax  ^  Ay.  The  figures  (2.1),  (2.2), 
and  (2.3)  visualize  the  local  error  for  the  (3-2),  (3-4),  and  (3-6)  FM  at  different  sections 
x  =  0.2, 0.4, 0.6, 0.8.  The  local  error  is  highest  at  extremum  points,  particularly  for  the 
second  order  method  which  is  globally  TVD.  However,  table  3  shows  that  the  global  order  of 
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accuracy  in  L1  and  L°°  norms  is  uniform.  Moreover,  table  4  shows  that  the  order  of  accuracy 
for  the  filtering  method  does  not  really  depend  on  the  value  of  the  parameter  a.  However, 


the  number  of  corrections  increases  as  cos  a  approaches  1  for  which  the  full  ENO  method  is 
performed. 


Scheme 


(3-2)FM 


(3-4)FM 


(3-6)FM 


Z^-norm  L°°-norm 


2.00 


1.45 


3.81 


3.35 


3.06 


3.11 


Table  3. 


cos  a 

#  of  corrections 

L1-norm 

L°°- norm 

CPU  time 

0.1  (3-4)FM 

84 

3.62 

3.01 

13.23 

0.9  (3-4)FM 

253 

3.61 

3.02 

13.45 

1.0  (3-4)RM 

1681 

3.91 

3.15 

36.41 

Table  4. 


From  this  table,  the  computational  cost  is  reduced  by  a  factor  of  almost  3  as  the  (3-4)FM 
is  used  (cos  a  <  0.9)  versus  the  full  ENO  (3-4)RF  method  (cos  a  =  1.).  /.Iso,  no  significant 
gain  in  the  order  of  accuracy  is  obtained  when  the  full  (3-4)RF  method  is  used. 


3.3  Example  III: 


In  this  example,  we  study  the  behavior  of  the  filtering  method  to  nonlinear  equations.  As 
simple  case,  we  consider  Burgers’  equation 


with  initial  condition 


u(x,  0)  =  cos  2IIx. 


We  again  take  periodic  boundary  conditions  on  the  interval  [0, 1].  The  solution  becomes 
discontinuous  at  time  t  =  ^  at  the  steady  location  x  =  0.25.  The  domain  ft  is  discretized 
and  the  grid  points  are  denoted  by  x*,  for  i  =  0,...,n,  and  let  n  =  40  in  the  numerical 
experiments.  The  exact  solution  is  approximated  by  using  Newton’s  method  whenever  the 
solution  is  smooth.  The  local  error  at  each  grid  points  is  plotted  for  the  (3-2),  (3-4),  and  (3- 
6)FM  in  the  figures  3.1,  3.2,  and  3.3.  Due  to  the  centered  differences  which  are  performed  in 
the  basic  scheme,  the  local  error  propagates  symmetrically  with  respect  to  inflection  points. 
Nevertheless,  the  table  5  shows  that  the  order  of  accuracy  of  the  filtering  method  is  still 
orm.  These  calculations  were  performed  at  time  t  =  0.1  with  a  CFL  coefficient  X  =  1.. 


uni 


Scheme 

Ll- norm  Z^-norm 

(3-2)FM 

1.95  1.40 

(3-4)FM 

3.74  3.69 

(3-6)FM 

3.50  3.47 

Table  5. 
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At  time  t  =  ^  the  shock  wave  appears.  The  shape  of  the  solution  is  shown  in  the  set 
of  figures  3.4,  3.5,  and  3.6.  No  spurious  oscillations  can  be  detected  from  tnese  plots.  The 
number  of  corrections  was  never  higher  than  8  whatever  the  number  of  grid  points.  Again 
no  visible  change  in  the  solution  occurs  when  different  values  of  a  are  considered.  In  this 
numerical  example,  cos  a  =  0.1  and  the  extremum  test  was  enforced. 

Using  Burgers’  equation  again,  we  want  to  test  whether  our  method  is  stable.  We  consider 
the  Riemann  problem  with  ui  =  +1,  and  uT  =  —  1  and  run  our  program  using  centered 
differences  only  on  a  40  point  grid  with  a  CFL  number  of  0.5  up  to  t  =  3.  The  figure  (3.7.1) 
visualizes  the  solution  at  this  time.  Large  oscillations  can  be  seen  up  to  the  boundaries. 
Taking  the  final  solution  of  the  previous  problem  as  initial  condition  and  using  the  (3-2)FM 
with  the  parameter  cos  a  =  0.99,  we  obtain  the  steady  solution  ui  =  +1,  ur  =  —1  (figure 
(3.7.2).  Moreover,  the  test  for  an  extrema  was  not  used.  This  shows  that  our  filtering  method 
acts  like  a  viscosity  method  in  regions  of  smoothness  even  as  the  parameter  cos  a  approaches 
1.  However,  the  number  of  corrections  occured  at  only  10  grid  points  on  the  average. 

Finally,  we  want  to  test  whether  our  method  works  for  nonconvex  fluxes  and  for  initial 
(  editions  solving  a  Riemann  problem  having  an  expansion  shock.  First  of  all,  we  considered 
Burgers’  equation  with  the  initial  condition  uj  =  —2  for  x  <  0.,  and  uT  =  2  for  x  >  0.  The 
centered  differences  of  the  basic  scheme  let  the  initial  expansion  shock  unchanged.  However, 
by  adding  a  small  perturbation  of  amplitude  e  =  10~3  to  the  initial  condition,  the  numerical 
solution  does  not  violate  the  entropy  condition  and  does  tend  for  long  time  to  the  stationary 
solution  of  the  problem,  i.e  u  =  0  (see  the  figure  3.7.3).  In  the  other  hand,  we  ran  our 
program  with  a  nonconvex  flux  f{u )  =  (u2~1Xu.3-4),  wjth  the  initial  condition  uj  =  2,  and 
uT  =  —2  in  each  side  of  x  =  0.  Again,  if  no  perturbation  is  added  to  this  initial  condition,  the 
centered  differences  let  the  solution  remain  unchanged.  The  results  with  a  small  initial  noise 
added  to  the  initial  condition  are  displayed  in  the  figure  3.7.4.  The  results  are  in  agreement 
with  those  presented  in  [6]. 

In  order  to  avoid  the  risk  of  developing  an  expansion  shock,  we  may  implement  another 
test  that  checks  whether  A^  <  0  <  Af ,  where  \f’R  are  the  eigenvalues  of  Vuf  in  each  side 
of  the  grid  point  x <  in  the  genuinely  nonlinear  fields  (for  non  convex  fields,  we  correct  at  all 
sonic  points).  If  this  test  is  satisfied  then  we  correct  the  numerical  solution  using  the  filter. 

3.4  Example  IV: 

In  this  example,  we  extend  example  III  to  two  space  dimensions  by  taking  the  two  dimen¬ 
sional  Burgers’  equation 
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with  the  same  initial  condition  as  in  example  II  and  1-periodic  boundary  conditions  in  both  x 
and  y  variables.  The  square  domain  is  discretized  as  in  example  II  with  m  =  40  and  n  =  30. 
The  shock  discontinuity  occurs  at  time  t  =  ^  at  the  steady  location  x  +  y  =  0.25.  We  study 
the  propagation  of  the  local  error  at  time  t  =  0.1  at  different  sections  x  =  0.2, 0.4, 0.6, 0.8. 
Again  the  error  is  symmetric  with  respect  to  inflection  points,  see  the  figures  (4.1),  (4.2),  and 
(4.3).  Moreover,  the  error  is  distributed  within  more  grid  points  as  the  order  of  the  space 
discretization  increases.  This  is  in  agreement  with  the  fact  that  the  width  of  the  stencil 
increases  with  the  order  of  accuracy.  Table  6  describes  the  order  of  accuracy  in  L 1  and 
L°°  norms  for  the  (3-2),  (3-4),  and  (3-6)  filtering  methods.  Table  7  discusses  the  order  of 
accuracy  in  these  norms  for  different  values  of  the  parameter  cos  a  for  the  (3,4)FM.  We  also 

M  method  versus  the  (3-4)LLF  ENO  method. 

Table  6. 

of  corrections  I^-norm 


211  3.24 

1681  3.33 

^b]e  ? 

The  CPU  time  is  again  reduced  by  a  factor  of  almost  3  when  the  filtering  method  is 
performed. 

At  time  ^  =  jn  ^e  shock  discontinuity  occurs  and  is  visualized  in  the  figure  4.4.  A  sharp 
transition  is  obtained  and  no  spurious  oscillations  can  be  seen.  Away  from  the  shock  the 
local  order  of  accuracy  is  preserved  showing  that  no  yopugation  of  error  starting  at  the 
shock  location  pollutes  the  smooth  part  of  the  numerical  solution.  Moreover,  the  formal 
orders  of  accuracy  are  approximately  the  same  as  those  given  in  the  table  6. 

3.5  Examples  V: 

In  these  examples,  we  study  the  smearing  of  contact  discontinuities  for  one  and  two  dimen¬ 
sional  problems.  To  illustrate  this  fact  and  study  the  behavior  of  the  filtering  method  in 
such  cases,  we  consider  the  linear  equations  given  in  examples  I  and  II  with  different  initial 
conditions: 

u(x,0)  =  cos(2IIx)  ±  0.5  If  0.25  <  x  or  x  >  0.75 
andu(x,0)  =  0  otherwise, 
u(x,y,0)  =  cos(2II(x  +y))  +  2  If  0.25  <  x,y  <  0.75 
andu(x,y,0)  =  0  otherwise, 


L°°-norm  CPU  time 
3T9  15.51 

3.25  18.32 

3.38  39.71 
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for  one  and  two  dimensional  problems,  respectively.  In  the  numerical  experiments  we  let 
n  =  40  and  study  the  numerical  solution  at  time  t  =  1.1  with  a  CFL  coefficient  A  =  0.5. 

For  the  one  dimensional  problem,  the  transition  is  plotted  in  the  set  of  figures  (5.1.1), 
(5.1.2),  and  (5.1.3)  for  the  (3-2),  (3-4),  and  (3-6)  FM,  respectively.  In  all  cases,  the  jump 
transition  is  localized  within  a  few  mesh  points.  However,  as  expected,  the  numerical  solution 
near  the  contact  discontinuities  is  better  approximated  for  the  highest  (sixth)  order  method. 
In  particular,  the  small  oscillations  that  appear  near  the  contact  discontinuities  for  the  second 
and  fourth  order  methods,  are  no  longer  there  for  the  sixth  order  method.  Also,  the  small 
oscillations  will  disappear  as  soon  as  the  value  of  the  parameter  cos  a  is  not  less  than  0.9. 
Moreover,  the  transition  from  the  upper  to  the  lower  parts  of  the  numerical  solution  is  better 
resolved  for  high  order  methods  (fourth  and  sixth  order),  whereas  the  transition  tends  to  be 
smeared  for  the  second  order  method. 

In  the  next  set  of  plots,  the  two  dimensional  problem  is  investigated.  The  results  are 
shown  in  the  figures  (5.2.1),  (5.3.1),  (5.4.1)  where  the  solution  wave  is  plotted  at  time 
t  =  1.1.  The  transition  from  the  zero  plateau  to  the  cosine  wave  is  smeared  within  three 
to  four  meshes  in  each  vertical  and  horizontal  sides  and  probably  more  at  each  of  the  four 
corners,  see  the  figures  (5.2.2),  (5.3.2),  and  (5.4.2).  Moreover,  the  maximum  error  occurs  at 
the  bottom  left  and  the  top  right  corners.  Also,  the  local  error  is  very  smooth  along  each 
sides  and  ’’discontinuous”  about  these  two  corners.  This  agrees  with  the  one  dimensional 
results  in  which  the  smearing  of  the  contact  discontinuities  happens  to  be  more  important 
in  the  upper  or  lower  part  of  the  cosine  wave. 

3.6  Examples  VI: 

The  last  examples  are  devoted  to  extend  our  filtering  method  to  systems  of  conservation 
laws.  We  consider  the  compressible  Euler  equations  for  gas  dynamics  introduced  during 
the  introduction.  The  Euler  equations  (  1)  and  (  4)  are  studied  with  these  sets  of  initial 
conditions: 

•  One  dimensional  Euler: 

—  We  consider  the  initial  condition  given  in  the  example  8  of  [7].  That  is,  we  take: 

p  =  3.857143,7  =  2.629369, p  =  10.3333333  when  x  <  —  4 
p  =  1  +  e  sin  5x,  7  =  0,p  =  1.  when  x  >  —4 

—  If  e  =  0.  we  get  a  pure  Mach  3  shock  moving  to  the  right.  Following  [7]  in 
example  8,  we  take  e  =  0.2  in  the  numerical  experiment. 
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•  Two  dimensional  Euler: 

—  We  consider  the  initial  condition  given  in  the  example  9  of  [7].  That  is,  we 
consider  a  Mach  8  shock  located  at  i  =  —  1  moving  to  the  right  into  the  state 
with 


Pr  =  1, 

Pr  =  1  > 

Or 

vx  = - sin  8r  cos  (xkrcos9r  +  ykr  sin  8r), 

Pr 

vv  =  —  cos  8t  cos  ( xkrcosOT  -f  ykT  sin  9r), 

Pr 

where  kr  —  ^  and  9r  =  j.  In  order  to  have  positive  pressure  during  the  calcula¬ 
tions,  we  used  a  parameter  cos  a  =  0.1  — >  a  ~  0.9y. 

The  results  for  the  one  dimensional  problem  are  shown  in  the  set  of  figures  (6.1.1),  (6.1.2), 
and  (6.1.3).  The  desired,  physical,  oscillations  near  the  shock  transition  which  are  parts  of 
the  exact  solution  are  particularly  visible  for  the  sixth  order  method.  The  second  order  and 
fourth  order  method  give  a  fairly  good  representation  of  the  expected  solution.  The  number 
of  corrections  was  about  25%  for  200  grid  points  and  a  CFL  coefficient  A  =  0.8. 

The  results  for  the  two  dimensional  shock- turbulence  problem  are  plotted  in  the  set  of 
figures  (6.2.1),  (6.2.2),  (6.3.1),  (6.3.2),  in  which  the  pressure  and  density  field  are  displayed. 
As  comparison,  similar  results  have  been  obtained  in  Using  this  value  of  the  parameter  cos  a, 
the  number  of  corrections  was  approximatively  40%  for  a  80  x  60  grid  with  a  CFL  coefficient 
of  0.5.  In  this  experiment,  the  (3-2)  and  (3-4)FM  have  been  used.  200  time  iterations  have 
been  necessary  to  reach  the  final  time  t  =  0.2.  70%  of  the  global  computational  time  was 
used  to  evaluate  the  ENO  interpolating  polynomials  during  the  filtering  step,  leading  to  a 
reduction  of  a  factor  of  2  of  the  computational  time  when  the  filtering  method  is  used. 

4  Final  Remarks  and  Conclusions 

The  main  conclusion  concerns  the  number  of  corrections  which  has  always  been  less  than 
10  to  25%  for  all  these  examples,  except  for  the  two  dimensional  shock-turbulence  problem 
in  which  a  very  complicated  structure  of  the  flow  appears.  Therefore,  the  computational 
cost  for  these  type  of  problems  is  reduced  by  a  factor  of  almost  3.  Moreover,  this  factor 
is  quite  significant  when  very  high  order  methods  are  implemented.  In  particular,  the  use 
of  sixth  order  method  for  the  one  dimensional  Euler  equations  leads  to  a  very  accurate 
approximate  solution.  An  important  remark  is  that  it  would  be  possible  in  the  near  future 
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to  implement  other  high  resolution  techniques  together  with  this  type  of  filtering  method, 
e.g  subcell  resolution  [4]. 

In  the  near  future,  this  type  of  method  will  be  implemented  in  for  different  type  of  prob¬ 
lems  involving  Hamilton-Jacobi  equations.  Moreover,  a  similar  approach  for  more  general 
domains  using  finite  element  triangulations  is  under  investigation.  So  far,  by  appropriate 
linear  combination  using  the  basis  functions  of  all  triangles  around  a  vertex,  it  has  been  pos¬ 
sible  to  construct  a  second  order  method  in  space  and  correct  the  oscillating  points  by  using 
a  first  order  Godunov  type  filter.  Higher  order  filtering  method  are  also  under  investigation. 
ENO  schemes  for  Hamilton-Jacobi  equations  in  Cartesian  coordinates  were  developed  in  [8] 
and  in  [9]. 
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Figure  Captions 

•  (2.1)  (3-2)FM:  x  =  0.2,  0.4, 0.6, 0.8  sections,  Local  Error  .10~2  For  2D  linear  Problem 
(Example  II). 

•  (2.2)  (3-4)FM:  x  =  0.2, 0.4, 0.6, 0.8  sections,  Local  Error  .10-5  For  2D  linear  Problem 
(Example  II). 

•  (2.3)  (3-6)FM:  x  =  0.2,  0.4, 0.6, 0.8  sections,  Local  Error  .IQ-5  For  2D  linear  Problem 
(Example  II). 

•  (3.1)  (3-2)FM:  Local  Error  .10~3  For  Burgers’  Equation  ( t  =  0.1). 

•  (3.2)  (3-4)FM:  Local  Error  .10-5  For  Burgers’  Equation  ( t  =  0.1). 

•  (3.3)  (3-6)FM:  Local  Error  .10~6  For  Burgers’  Equation  ( t  =  0.1). 

•  (3.4)  (3-2)FM:  Shock  Transition  for  Burgers’  Equation  at  time  t  = 

•  (3.5)  (3-4)FM:  Shock  Transition  for  Burgers’  Equation  at  time  t  =  jTi- 

•  (3.6)  (3-6)FM:  Shock  Transition  for  Burgers’  Equation  at  time  t  = 

•  (3.7.1)  (3-2)CD  (Centered  Differences):  Solution  Wave  of  Burgers’  Equation  at  time 
t  =  1.  with  +1,-1  Initial  Condition. 

•  (3.7.2)  (3-2)FM:  Solution  Wave  of  Burgers’  Equation  at  time  t  =  1.  with  Initial  Con¬ 
dition  of  Figure  (3.7.1). 

•  (3.7.3)  (3-4)FM:  Solution  Wave  of  Burgers’  Equation  at  time  t  —  l.,2.,3.  with  Initial 
Condition  tq  =  —  2,ur  =  2  plus  ±10“3  noise. 

•  (3.7.4)  (3-4)FM:  Solution  Wave  of  nonconvex  nonlinear  Equation  ut  +  (u3~xKu  *  =  q 

at  time  t  =  0.2, 0.4, 0.8  with  Initial  Condition  tq  =  2,ur  =  —2  plus  ±10-3  noise. 

•  (4.1)  (3-2)FM:  Local  Error  For  2D  Burgers’  Equation  (*  =  0.1). 

•  (4.2)  (3-4)FM:  Local  Error  For  2D  Burgers’  Equation  ( t  =  0.1). 

•  (4.3)  (3-6)FM:  Local  Error  For  2D  Burgers’  Equation  ( t  =  0.1). 

•  (4.4)  (3-2)FM:  Shock  Transition  For  2D  Burgers’  Equation  at  time  t  — 

•  (5.1.1)  (3-2)FM:  Linear  ID  Equation  -  Contact  Discontinuity  at  time  t  =  1.1. 
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•  (5.1.2)  (3-4)FM:  Linear  ID  Equation  -  Contact  Discontinuity  at  time  t  =  1.1. 

•  (5.1.3)  (3-6)FM:  Linear  ID  Equation  -  Contact  Discontinuity  at  time  t  =  1.1. 

•  (5.2.1)  (3-2)FM:  Linear  2D  Equation  -  Contact  Discontinuity  at  time  t  =  1.1. 

•  (5.2.2)  (3-2)FM:  Linear  2D  Equation  -  Error  Plot  and  Contour  -  Contact  Discontinuity 

at  time  t  =  1.1. 

•  (5.3.1)  (3-4)FM:  Linear  2D  Equation  -  Contact  Discontinuity  at  time  t  =  1.1. 

•  (5.3.2)  (3-4)FM:  Linear  2D  Equation  -  Error  Plot  and  Contour  -  Contact  Discontinuity 
at  time  t  =  1.1. 

•  (5.4.1)  (3-6)FM:  Linear  2D  Equation  -  Contact  Discontinuity  at  time  t  =  1.1. 

•  (5.4.2)  (3-6)FM:  Linear  2D  Equation  -  Error  Plot  and  Contour  -  Contact  Discontinuity 
at  time  t  =  1.1. 

•  (6.1.1)  (3-2)FM:  ID  Euler  Equations  e  =  0.2. 

•  (6.1.2)  (3-4)FM:  ID  Euler  Equations  c  =  0.2. 

•  (6.1.3)  (3-6)FM:  ID  Euler  Equations  e  =  0.2. 

•  (6.2.1)  (3-2)FM:  2D  Euler  Equations  -  Shock-Turbulence  Interaction  -  Pressure  field 
t  =  0.2. 

•  (6.2.2)  (3-2)FM:  2D  Euler  Equations  -  Shock-Turbulence  Interaction  -  Density  field 

t  =  0.2. 

•  (6.3.1)  (3-4)FM:  2D  Euler  Equations  -  Shock  Turbulence  Interaction  -  Pressure  field 
t  =  0.2. 

•  (6.3.2)  (3-4)FM:  2D  Euler  Equations  -  Shock  Turbulence  Interaction  -  Density  field 
t  =  0.2. 
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- 1=0. 1  -  -  Local  Error  at  time  1=0. 1  - 


-  (3,6)  CD  scheme  +  uniform  EN O-Roe -6  filter  -  -  figure  2.3  - 


xlO3  -  (3,2)  filter  scheme  -  Error  LI  =  1.95  ,  Error  Linf  =  1.40  - 


-  Pointwise  Error  ■  Figure  3. 1  - 
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-  t=l/(2*pi)  -  -t=l/(2*pi)  - 
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-  t=l/(2*pi)  - 


•  (3,6)  filter  scheme  ,  Burger  equation  with  cosine  wave  - 


-  Solution  wave  -  Figure  3.7.2  - 
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-  (3,6)  fill a  acheme ,  Linar  equation  .Contact  discontinuity  - 


Line  »r  Equation  Ui+(U)x+(U)y=0  .Contact  discontinuity 
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